C
C     the incident field is expanded in as0, bs0
C     PW case and DIP cases are both workable
C     
C     The main function of this subroutine is solving 
C     the interaction matrix with bicgstab
C
c  ------------------------------------------------------------
c  begins iteration process to solve the interaction equations
c  for partial interacting scatteing coefficients
c  factor1,factor2=0: BI-CGSTAB [see H.A. van der Vorst, SIAM,
c  J. Sci. Stat. Comput. 13, 631 (1992)]
c  factor1=factor2=1: order-of-scattering [Fuller and Kattawar,
c  Opt. Lett. 13, 90 & 1063 (1988)]
c  0<factor1,factor2<1: an iteration scheme [Xu, Appl. Opt. 34,
c  4573 (1995)]
c  When the number of iterations exceeds the maximum MXINT, it
c  switches to use BI-CGSTAB.
c  In the iteration solution of interacting equations, a vector
c  of wave expansion referred to a sphere-center needs to be
c  translated to other sphere-centers, which demands the storage
c  and evaluation of a large number of vector translation
c  coefficients. This code uses Mackowski's three-step technique
c  (rotation-translation-rotation method) [see Mackowski, Proc.
c  R. Soc. Lond. A 433, 599 (1991)] for the translation process,
c  which significantly reduces both computer memory usage and
c  computing time requirements by decomposition of the vector
c  translation coefficients into rotational and axial
c  translational parts. By using Mackowski's three-step method,
c  the calculations involve only translation along z-axis.
c  -------------------------------------------------------------
      subroutine Mie_interaction(nL, as,bs, aMie,bMie, p0,q0, r0, nmax,
     $     uvmax, ek, drot, atr0, btr0, factor1, factor2, MXINT,fint
     $     ,small)
      implicit double precision(a-h, o-z)
      include 'gmm01f.par'
      parameter (nmp=np*(np+2),nmp0=(np+1)*(np+4)/2)
      
      integer ind(nLp), nmax(nLp), uvmax(nLp)

      double precision c0i(nLp), c1i(nLp), r0(6,nLp)

      complex*16 A,B,cmz,Aj,Bj,A2,B2,Aj2,Bj2,A0,B0,ephi,ci,cin,
     +     p0(nLp,nmp), q0(nLp,nmp), aMie(nLp,np),bMie(nLp,np),
     +	   as(nLp,nmp),bs(nLp,nmp),  as0(nLp,nmp),bs0(nLp,nmp),
     +     asc(nLp,nmp),bsc(nLp,nmp),as1(nLp,nmp),bs1(nLp,nmp),
     +     ast(nLp,nmp),bst(nLp,nmp),asp(nLp,nmp),bsp(nLp,nmp),
     +     asv(nLp,nmp),bsv(nLp,nmp)

C     only one sphere, do not need to solve interaction problem
      if(nL.eq.1) goto 490

      factor = factor1

      do i = 1,nL
         ind(i) = 0
         c0i(i) = 0.d0
         do n = 1,nmax(i)
            imn = n*n + n + 1
            c0i(i) = c0i(i) + p0(i,imn)*dconjg(p0(i,imn))
            c0i(i) = c0i(i) + q0(i,imn)*dconjg(q0(i,imn))
            c0i(i) = c0i(i) + p0(i,imn-2)*dconjg(p0(i,imn-2))
            c0i(i) = c0i(i) + q0(i,imn-2)*dconjg(q0(i,imn-2))
         enddo
      enddo

      niter = 1
      if(factor1.lt.0.001d0.or.factor2.lt.0.001d0) goto 61
      write(6,*) 'Starting iteration solution process'
      do i = 1,nL
         do imn = 1,uvmax(i)
            as0(i,imn) = p0(i,imn)
            bs0(i,imn) = q0(i,imn)
         enddo
      enddo

 60   call trans(nL,r0,nmax,uvmax,fint,atr0,btr0,ek,drot,as0,bs0,as1
     $     ,bs1,ind)

      do i = 1,nL
         if(ind(i).gt.0) goto 601
         c1i(i)=0.d0
         do imn=1,uvmax(i)
            n=dsqrt(dble(imn))
            as0(i,imn) = p0(i,imn) - aMie(i,n)*as1(i,imn)
            bs0(i,imn) = q0(i,imn) - bMie(i,n)*bs1(i,imn)
            A = as0(i,imn) - as(i,imn)
            B = bs0(i,imn) - bs(i,imn)
            c1i(i) = c1i(i) + A*dconjg(A)
            c1i(i) = c1i(i) + B*dconjg(B)
            as0(i,imn) = as(i,imn) + factor*A
            bs0(i,imn) = bs(i,imn) + factor*B
            as(i,imn) = as0(i,imn)
            bs(i,imn) = bs0(i,imn)
         enddo
 601     continue
      enddo

C     check convergence
      cext0 = 0.d0
      cext1 = 0.d0
      do i = 1,nL
         if(ind(i).gt.0) goto 602
         cext0 = cext0 + c0i(i)
         cext1 = cext1 + c1i(i)
         temp = c1i(i)/c0i(i)
         if(temp.lt.small) ind(i)=1
 602     continue
      enddo
C     when convergence, finish iteration by go to 490
      temp = cext1/cext0
      if(temp.lt.small) goto 490
      write(6,'(a11,i4,2x,e15.7)') 'iteration #', niter, temp

C     when reaches the maximun iteration switch to bicg
      if(niter.gt.MXINT) then
         write(6,*) '  *** Maximum iterations exceeded ***'
         write(6,*) '  *** Switched to Bi-CGSTAB method***'
         do i = 1,nL
            ind(i) = 0
            do imn = 1,uvmax(i)
               as(i,imn) = p0(i,imn)
               bs(i,imn) = q0(i,imn)
            enddo
         enddo
         niter = 1
         goto 61
      endif
C     continue another iteration by go to 60
      niter = niter + 1
      goto 60

C     entry of bicg
 61   write(6,*) 'Starting Bi-CGSTAB solution process'
      call trans(nL,r0,nmax,uvmax,fint,atr0,btr0,ek,drot,as,bs,as1,bs1
     $     ,ind)
      do i=1,nL
         c1i(i)=0.d0
         do imn=1,uvmax(i)
            n = dsqrt(dble(imn))
            as1(i,imn) = -aMie(i,n)*as1(i,imn)
            bs1(i,imn) = -bMie(i,n)*bs1(i,imn)
            c1i(i) = c1i(i) + as1(i,imn)*dconjg(as1(i,imn))
            c1i(i) = c1i(i) + bs1(i,imn)*dconjg(bs1(i,imn))
         enddo
      enddo
      temp = 0.d0
      do i=1,nL
         cext0 = c1i(i)/c0i(i)
         if(cext0.lt.small) ind(i)=1
         if(cext0.gt.temp) temp = cext0
      enddo
      if(temp.lt.small) goto 490
      A0 = 0.d0
      do i=1,nL
         if(ind(i).gt.0) goto 611
         do imn=1,uvmax(i)
            asp(i,imn) = as1(i,imn)
            bsp(i,imn) = bs1(i,imn)
            as0(i,imn) = as1(i,imn)
            bs0(i,imn) = bs1(i,imn)
            A0 = A0 + as1(i,imn)*as1(i,imn)
            A0 = A0 + bs1(i,imn)*bs1(i,imn)
         enddo
 611     continue
      enddo
C
 62   call trans(nL,r0,nmax,uvmax,fint,atr0,btr0,ek,
     +     drot,asp,bsp,asv,bsv,ind)
      do i=1,nL
         if(ind(i).gt.0) goto 621
         do imn=1,uvmax(i)
            n=dsqrt(dble(imn))
            asv(i,imn)=aMie(i,n)*asv(i,imn)+asp(i,imn)
            bsv(i,imn)=bMie(i,n)*bsv(i,imn)+bsp(i,imn)
         enddo
 621     continue
      enddo
      A=0.d0
      do i=1,nL
         if(ind(i).gt.0) goto 622
         do imn=1,uvmax(i)
            A=A+asv(i,imn)*as1(i,imn)
            A=A+bsv(i,imn)*bs1(i,imn)
         enddo
 622     continue
      enddo
      Aj = A0/A
      do i=1,nL
         if(ind(i).gt.0) goto 623
         do imn=1,uvmax(i)
            asc(i,imn) = as0(i,imn)-Aj*asv(i,imn)
            bsc(i,imn) = bs0(i,imn)-Aj*bsv(i,imn)
         enddo
 623     continue
      enddo
      call trans(nL,r0,nmax,uvmax,fint,atr0,btr0,ek,drot,asc,bsc,ast
     $     ,bst,ind)
      do i = 1,nL
         if(ind(i).gt.0) goto 624
         do imn=1,uvmax(i)
            n=dsqrt(dble(imn))
            ast(i,imn) = aMie(i,n)*ast(i,imn) + asc(i,imn)
            bst(i,imn) = bMie(i,n)*bst(i,imn) + bsc(i,imn)
         enddo
 624     continue
      enddo
      A2=0.d0
      B2=0.d0
      do i=1,nL
         if(ind(i).gt.0) goto 625
         do imn=1,uvmax(i)
            A2=A2+ast(i,imn)*asc(i,imn)
            A2=A2+bst(i,imn)*bsc(i,imn)
            B2=B2+ast(i,imn)*ast(i,imn)
            B2=B2+bst(i,imn)*bst(i,imn)
         enddo
 625     continue
      enddo
      Bj = A2/B2
      do i=1,nL
         if(ind(i).gt.0) goto 626
         c1i(i)=0.d0
         do imn=1,uvmax(i)
            Aj2=Aj*asp(i,imn)+Bj*asc(i,imn)
            Bj2=Aj*bsp(i,imn)+Bj*bsc(i,imn)
            c1i(i)=c1i(i)+Aj2*dconjg(Aj2)
            c1i(i)=c1i(i)+Bj2*dconjg(Bj2)
            as(i,imn) = as(i,imn)+Aj2
            bs(i,imn) = bs(i,imn)+Bj2
         enddo
 626     continue
      enddo

c     bicg迭代收敛测试
      cext0 = 0.d0
      cext1 = 0.d0
      do i=1,nL
         if(ind(i).gt.0) goto 627
         cext0 = cext0 + c0i(i)
         cext1 = cext1 + c1i(i)
 627     continue
      enddo
      temp = cext1/cext0
      if(temp.lt.small) goto 490
      if(niter.gt.MXINT) then
         write(6,*) 'Caution:'
         write(6,*) '  *** Maximum iterations exceeded ***'
         goto 490
      endif

      write(6,'(a11,i4,2x,e15.7)') 'iteration #',niter,temp

      do i=1,nL
         if(ind(i).gt.0) goto 628
         do imn=1,uvmax(i)
            as0(i,imn) = asc(i,imn) - Bj*ast(i,imn)
            bs0(i,imn) = bsc(i,imn) - Bj*bst(i,imn)
         enddo
 628     continue
      enddo
      A2 = 0.d0
      do i=1,nL
         ast(i,1) = 0.d0
         if(ind(i).gt.0) goto 629
         do imn=1,uvmax(i)
            ast(i,1) = ast(i,1) + as0(i,imn)*as1(i,imn)
            ast(i,1) = ast(i,1) + bs0(i,imn)*bs1(i,imn)
         enddo
         A2 = A2 + ast(i,1)
 629     continue
      enddo
      B0 = A2/A0
      B0 = B0*Aj/Bj
      do i=1,nL
         if(ind(i).gt.0) goto 630
         do imn=1,uvmax(i)
            asp(i,imn) = as0(i,imn) + B0*(asp(i,imn)-Bj*asv(i,imn))
            bsp(i,imn) = bs0(i,imn) + B0*(bsp(i,imn)-Bj*bsv(i,imn))
         enddo
 630     continue
      enddo
      A0=0.d0
      do i=1,nL
         if(ind(i).gt.0) goto 631
         cext0 = c1i(i)/c0i(i)
         if(cext0.lt.small) then
            ind(i)=1
            goto 631
         endif
         A0 = A0 + ast(i,1)
 631         continue
      enddo
C     进入下一次bicg迭代
      niter = niter+1
      goto 62

c ---------------------------------------------------------------
c Solution of scattering coefficients finished. Label 490 is where
c the solver will jump on convergence or error.
c ---------------------------------------------------------------
 490  do i = 1,nL
         ind(i) = 0
      enddo

c ---------------------------------------------------------------
c Writing scattering coefficients.
c ---------------------------------------------------------------
c	open(17,file='scacof.dat',status='unknown')
c	do 282 i=1,nL
c	   write(17,'(a2,i3)') 'i=',i
c	   do 281 j=1,uvmax(i)
c	      n=dsqrt(dble(j))
c	      m=j-n*n-n
c	write(17,'(2i5,4e17.8)') m,n,dble(as(i,j)),dimag(as(i,j)),
c     +                               dble(bs(i,j)),dimag(bs(i,j))
c 281      continue
c 282   continue
c	close(17)

      write(6,*) 'Scattering coefficients done.'

      return
      end
